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The accuracy and stability of the Caltech-Cornell pseudospectral code is evaluated using the 
Kidder, Scheel, and Teukolsky (KST) representation of the Einstein evolution equations. The basic 
"Mexico City Tests" widely adopted by the numerical relativity community are adapted here for 
codes based on spectral methods. Exponential convergence of the spectral code is established, 
apparently limited only by numerical roundoff error or by truncation error in the time integration. 
A general expression for the growth of errors due to finite machine precision is derived, and it is 
shown that this limit is achieved here for the linear plane- wave test. 
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PACS numbers: 04.25.Dm, 02.70.Hm, 02.60.Cb 
I. INTRODUCTION 

A number of groups have now developed numerical rel- 
ativity codes sophisticated enough to evolve binary black- 
hole spacetimes 0, 0, i, i, i, 6] . The gravitational wave- 
forms predicted by these evolutions will play an impor- 
tant role in detecting and interpreting the physical prop- 
erties of the sources of these waves soon to be detected 
(we presume) by the community of gravitational wave 
observers (e.g., LIGO, etc.). Therefore, such codes must 
be capable of performing stable and accurate simulations 
of very nonlinear and dynamical spacetimes. 

Several years ago a large subset of the numerical rel- 
ativity community, the "Apples with Apples" collabora- 
tion [31 , proposed a scries of basic code tests designed to 
verify the accuracy, stability, robustness and efficiency of 
any code designed to find fully three-dimensional solu- 
tions to the Einstein evolution equations. These tests — 
often referred to as the "Mexico City Tests" because they 
were first formulated during a conference in Mexico City 
in May 2002 — were designed to be analogous to the stan- 
dard suite of tests used by the numerical hydrodynam- 
ics community (e.g., tests to reproduce Sedov explosions. 
Sod shock tubes, blast waves, etc.) to commission new 
hydrodynamics codes. The Mexico City tests were de- 
signed to be applicable to any formulation of Einstein's 
equations solved with any numerical method. All tests 
proposed so far concern bulk properties of the formu- 
lation and numerical method, and so all of the evolu- 
tions are carried out on a numerical grid with three-torus 
topology; no boundary conditions are needed (or tested). 
There arc four basic tests, some of them in a number 
of variations: (a) the evolution of certain "random" ini- 
tial data; (b) the evolution of small-amplitude "linear" 
plane-wave initial data; (c) the evolution of a nonlinear 
gauge-wave representation of fiat spacetime; and (d) the 
evolution of initial data for a very dynamic and nonlinear 
Gowdy cosmological model. 

The Mexico City tests have now been applied to a num- 
ber of different numerical relativity codes that use differ- 
ent formulations of the Einstein equations 0, [l] ■ But all 



of the codes tested so far use finite difference numerical 
methods. In this paper we report the results of applying 
these tests to the code being developed by the collabo- 
ration between the Caltech and Cornell numerical rela- 
tivity groups. We use a first-order symmetric-hyperbolic 
formulation of the equations developed by Kidder, Scheel 
and Teukolsky Q (sometimes referred to as the KST for- 
mulation) and we solve the equations using pseudospec- 
tral numerical methods. The results reported here dif- 
fer therefore from all previously tested cases both in the 
formulation of the Einstein equations and the numerical 
methods used to solve them. 

This paper is organized as follows: In Sec. [ll] we re- 
view the KST formulation of the equations and the pseu- 
dospectral numerical methods that we use to solve them. 
The remaining sections present the results of the various 
Mexico City tests, adapted somewhat to provide more 
challenging tests of a code based on spectral methods. 
In Sec, mil we show that our code is stable when evolving 
small random perturbations of fiat spacetime. In Sec. IIVI 
we report the results of the small-amplitude plane-wave 
test. We demonstrate the convergence rates for differ- 
ent spatial resolutions and different time-step algorithms. 
We also derive an equation for the error introduced by 
finite machine precision, and show that it limits the con- 
vergence of our evolutions for small spatial and tempo- 
ral resolutions. In Sec. |V] we show the stability of our 
evolution code for nonlinear gauge waves. In this case, 
nonlinear terms give rise to an instability that is drasti- 
cally reduced by suitably filtering the components of the 
spectral expansion. Section IVTl shows the performance of 
our code for evolutions of the highly dynamical Gowdy 
spacetime, in which the exact analytical expressions for 
the components of the fields grow exponentially in time. 
Finally, we discuss and summarize our various results in 

Sec. inn 



II. SOLUTION METHOD 

In this section we describe the formulation of the Ein- 
stein equations and the pseudospectral numerical solu- 
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tion method that we test. The Mexico City tests were 
designed with finite-difFerence methods in mind and were 
originally applied to formulations of the Einstein equa- 
tions that are second-order in space and first-order in 
time. Both our numerical methods and our represen- 
tation of the Einstein equations differ significantly from 
those in Ref . 0] , so appropriate modifications to the Mex- 
ico City test suite (for example, the number of grid points 
used or the constraint quantities observed) are needed. 
These modifications are also described in this section. 



A. KST Formulation 

The KST system is a first-order symmetric hyper- 
bolic generalization of York's representation of the ADM 
equations The dynamical variables of this system 

are the three- metric gij , the extrinsic curvature Kij , and 
a new variable Dkij that is initially set equal to OkSij/'i. 
This last variable allows the system to be put into first- 
order form. Its introduction results in two additional 
constraints: 
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The KST evolution equations are obtained from the 
ADM equations [l^l by adding constant multiples of the 
various constraints to the evolution equations and by re- 
placing the lapse with a lapse density function. These 
changes do not affect the physical solutions of the system, 
but they do modify the unphysical constraint-violating 
solutions. The added constraint terms are proportional 
to constant parameters {71,72,73,74}, which are cho- 



sen to make the system symmetric hyperbolic The 
principal parts of the KST evolution equations, then, are 
given by: 
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Here, the symbol ~ indicates that terms algebraic in the 
fields (that is, nonprincipal terms) are not shown explic- 
itly. The lapse function N is taken to be 



N 



(6) 



and both the lapse density function Q and the shift N'^ 
are assumed to be specified functions of the coordinates. 



rather than independent dynamical fields. Since each 
of the Mexico City tests involves reproducing either a 
known analytic solution of the Einstein equations or a 
small perturbation about a known solution, for all tests 
reported here we set the lapse density Q and the shift TV' 
from the appropriate analytic solution. We choose one 
set of the KST parameters for all the tests here: 70 = 0.5; 
71 = -0.21232; 72 = -0.00787402; 73 = -I.6I994; 74 = 
—0.69885. These values were chosen because they make 
the KST system symmetric hyperbolic and coincide with 
a set preferred by Owen [l^l in his extension of the KST 
system. 

To evaluate errors it is useful to look at constraint 
quantities. As mentioned above, the KST system has ad- 
ditional constraints, Eqs. ([T]) and ([2]), besides the usual 
Hamiltonian constraint C and momentum constraint Ci. 
To ensure that we are satisfying all the constraints, we 
monitor a single quantity C that is zero if and only if all 
of the constraints vanish: 



C 



(7) 



where an object is squared using the evolved spatial met- 
ric: for example, (C^)^ = g^^CiCj. 

Likewise, when evaluating differences from analytically 
known solutions, we define an overall error quantity that 
includes the errors in all evolved variables gij, Kij, and 
Dk^J. Taking 5g,j = - and similarly 

for other fundamental fields, this overall error quantity 
is given by 



6U = ^ {5g,jf + {5K,jf + {6Dk^Jf . (8) 

Notice that 5 U vanishes if and only if all evolved variables 
match the known solution. 

For all error quantities Q we display L2 norms: 



2|l2^ / Q'^/\9\dH 



(9) 



where Vol = J yJ\g\(Px is the volume of the domain. 
These norms are computed after each time step over the 
current t = const, hypersurface. We refer to ||C||2 as the 
constraint energy, and ||(5W||2 as the error energy. 

The error quantities ||(5W||2 and j|Cj|2 scale with the 
absolute magnitude of the fundamental fields and their 
derivatives, so it can be difficult to judge the significance 
of these error measures without knowing the overall scale 
of the variables in the problem. For this reason, we some- 
times plot the normalized error energy ||(5W||2/||Z^||2 and 
the normalized constraint energy ||C||2/||9Z-/||2; where the 
normalization factors are defined by 

U ^ ^ {g,^f + (K^^f + {D,i,f , (10) 
dU = \/ {d.,gjkf + + {d,D,kif ■ (H) 
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Note that \\5U\\2/\\U\\2 and ||C||2/||5W||2 become of order 
unity when errors completely dominate the numerical so- 
lution. We display normalized error quantities only for 
tests involving the Gowdy spacetimes fSection IVip . in 
which the fundamental variables vary exponentially in 
time. All other tests presented here involve perturba- 
tions of Minkowski spacctime, in which case the quantity 
||9Z-/||2 is of order the size of the perturbation and is 
therefore inappropriate to use as a normalization factor. 
However, for perturbations of Minkowski spacetime, the 
overall scale is of order unity so it suffices to display the 
unnormalized quantities ||(5ZY||2 and ||C||2. 

B. Pseudospectral Methods 

All of our numerical computations are carried out using 
pseudospectral methods; this is the first time the Mexico 
City tests have been applied to a pseudospectral code. A 
brief outline of our method is as follows: Given a system 
of partial differential equations 

9tw(x, t) = T[u{^, t), a,u(x, t)] , (12) 

where u is a collection of dynamical fields, the solution 
m(x, t) is expressed as a time-dependent linear combina- 
tion of TV spatial basis functions (j>k{^)- 

N-l 

u{x,t) = {tfe(i)(?!)fc(x) . (13) 

k=0 

Associated with the basis functions is a set of Nc col- 
location points Xi. Given spectral coefficients Uk{t), the 
function values at the collocation points u(yii,t) are com- 
puted using Eq. ([T^ . Conversely, the spectral coefficients 
are obtained by the inverse transform 

Uk{t) ^ ^ Wiu{xi,t)(j)k{xt) , (14) 

i=0 

where Wi are weights specific to the choice of basis func- 
tions and collocation points. Thus it is straightforward 
to transform between the spectral coefficients Ukit) and 
the function values at the collocation points u{xi,t). 

To solve the differential equations, we evaluate spatial 
derivatives analytically using the known derivatives of 
the basis functions, 

AT-l 

diu{x,t) = ^ Uk{t)di(j)k{x) , (15) 

k=0 

and we evaluate nonlinear terms using the values of 
u(xj,t) at the collocation points. Thus we can write the 
partial differential equation, Eq. p^ . as a set of ordi- 
nary differential equations for the function values at the 
collocation points, 

dtu{xi,t) = g^[u{-Kj,t)] , (16) 



where Qi depends on w(xj, t) for all j. We then integrate 
this system of ordinary differential equations in time, us- 
ing (for example) a fourth-order Runge-Kutta algorithm. 

Because the tests discussed here are periodic in all spa- 
tial dimensions, we use Fourier basis functions. If we 
choose a computational domain extending from —1/2 to 
1/2 in each of the x-, y-, and z-dircctions, then each vari- 
able u is decomposed as 

N^-l Ny-l 

'u.{x,y,z) = X! X! X! o-kini<l>k{x)(t)i{y)(t)„i{z), 

k=0 1=0 m=0 

(17) 

where 

r 1 /c = ; 

0fe(x) = < sin[7ra;(fc + 1)] k>0 (fc odd) ; (18) 
I cos(7rxfc) fc > (k even) . 

For smooth solutions, the spectral approximation 
Eq. converges exponentially (error ~ e"^''^ for some 
A > which depends on the solution). This is much 
faster than the polynomial convergence (error ~ l/N^) 
obtained using pth-order finite-differencing. As a result, 
we run our tests at coarser resolutions than those recom- 
mended in Ref. 01 for finite-difference codes — typically 
we use Ni = 9, 15, 21, 27, and 33 collocation points in 
the relevant directions. From Eqs. ([T7)) and ^TE\\ we see 
that if we choose N^, Ny, or to be an even integer, 
the highest frequency component in our expansion will 
have a sine term but no matching cosine term. Con- 
sequently the spatial derivative of this highest frequency 
component will not be represented by our basis functions, 
causing a numerical instability. Therefore we choose N^, 
Ny, and Nz to be odd. 

Because spectral methods so greatly reduce spatial 
discretization errors, time-stepping error is often domi- 
nant. In order to make the time stepping and the spatial 
discretization errors comparable in these tests, we use 
fourth-order Runge-Kutta ODE integration. The time- 
step sizes are chosen in an effort to use step sizes com- 
parable to those used to test finite difference methods 
in Ref. 01, while also ensuring that time-step errors do 
not dominate over our spatial truncation errors. We use 
At = Ax/20 in the first test, and At = Ax/AO in all 
others, except where explicitly noted. Here, Ax is the 
minimum distance between collocation points. 

III. RANDOM INITIAL DATA ON FLAT SPACE 

Perhaps the simplest test of a numerical relativity code 
is evolving standard Minkowski spacetime on a three- 
torus, T^. However, this test is too simple because all 
fundamental fields are spatially constant and most are 
identically zero, and hence most numerical methods will 
reproduce the correct solution exactly. This test can be 
made more discriminating by adding a small amount of 
random noise to the initial data; the noise is intended to 
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FIG. 1: Constraints for Minkowski with Random 
Noise. Higher resolutions are expected to have larger con- 
straints because more closely spaced points result in larger 
derivatives. The constraints do not grow in time. 



simulate the effect of finite numerical precision. A differ- 
ent random number is added to each component of each 
evolved variable, at each point in the domain. These 
random numbers are chosen to lie between —10^^" and 
10"^*^ so that the system remains in the linear regime. 
If these small perturbations to a simple spacetime grow 
unstably, it is likely that the inevitable errors (e.g., dis- 
cretization error or even numerical roundoff error) that 
arise in any more complicated simulation will also grow 
unstably. For this test we vary the resolution in the x- 
dimension, and we fix the resolution to three collocation 
points in each of the y- and z-dimensions. 

If the perturbations in the fields are chosen to be of size 
e independent of resolution, then the perturbation in the 
nth spatial derivatives of these fields will be ~ e{Ax)~", 
where Aa; is some measure of the distance between neigh- 
boring points. This means that error quantities involving 
derivatives (such as constraints) will be larger for finer 
resolutions.^ This behavior is seen in the plot of the con- 
straint energy in Fig. [1] 

The purpose of this test is to establish that small con- 



^ The Mexico City tests collaboration Q intended their Hamil- 
tonian constraint errors to be resolution-independent, so they 
chose the size of the perturbation e to be resolution- dependent, 
e ~ (Ax)^, which is the appropriate scaling for the second-order- 
in-space formulations of Einstein's equations they use. However, 
for the first-order-in-spacc formulation we use, the Hamiltonian 
constraint is computed using first derivatives of D^ij rather than 
second derivatives of gij, so the constraint will vary as (Aa;)"'^. 
Note also that the e ~ (Ax)-^ scaling docs not make the momen- 
tum constraint independent of resolution, as it depends on first 
derivatives of the fields. Wo simply choose e to be independent 
of resolution. 
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FIG. 2: Error Energy for Minkowski with Random 
Noise. The linear increase in time is due to a nonzero average 
in the random noise added to Kij. This average approaches 
zero as resolution is increased, since there are more points over 
which to average. The flat line shows the evolution when the 
average value of Kij is set to zero in the initial data. 



straint violations around flat space do not grow, and the 
KST system clearly passes this test. Whether or not 
constraint violations decay will depend on the evolution 
system and the numerical method. For example, arti- 
ficial dissipation in the numerical method might cause 
all variations to decay, including constraint violations. 
Furthermore, if the evolution system contains constraint 
damping in some form, then the constraints should de- 
cay. Indeed, Owen has extended the KST system to in- 
clude constraint damping [T^ : running the same test, he 
observes exponential decay in the constraint quantities. 
The flat constraint violations observed in Fig. [1] indicate 
that the KST system with our parameter choice does not 
damp constraints and that the spectral method has in- 
significant artificial dissipation. 

In Fig. [2] we see a linear growth of the error energy 
\\SU\\ for this test. We find that the growth is caused 
solely by contributions from the metric g^ ; the average 
values of Kij and Dkij remain constant in time. We can 
understand this as follows: The average value of Kij is 
determined by the random initial data and will in gen- 
eral be nonzero. The time derivative of Kij , to first order 
in the amplitude of perturbations around flat space, in- 
volves only spatial derivatives of Dkij ■ These derivatives 
have zero average (up to roundoff errors ~ 10~^^), be- 
cause the constant term in the Fourier expansion Eq. (jlSp 
is removed by differentiation, and therefore the average 
of Kij will be constant in time. The time derivative of 
gij involves a term proportional to Kij . Because the av- 
erage of Kij is constant in time and nonzero, the value 
of gij will therefore drift linearly in time. The average of 
Kij is smaller for higher resolutions — because one aver- 



5 



ages over more random numbers — which means that the 
growth rate of gij should decrease with increasing reso- 
lution. Indeed, this is what we observe in Fig. [21 

We can verify that the nonzero average of Kij is the 
only source of growth in gij by manually removing the 
average value of Kij . We expect this will leave the norms 
of the components of gij approximately constant in time. 
This is accomplished by setting the fc = spectral coeffi- 
cients of all components of Kij to zero in the initial data, 
after all the random numbers have been added. The flat 
line in Fig.[2]shows the result, indicating that the average 
offset in Kij is the only source of growth in the evolved 
variables of the KST system for this test. 



IV. LINEAR PLANE WAVE 

If the ultimate goal of simulating binary black hole 
mergers is to predict the gravitational-radiation wave- 
forms for observations, an evolution system must at least 
be capable of evolving a simple linear plane wave through 
flat spacetime. The form suggested for the Mexico City 
tests in Ref. is 

ds^ = -dt^ + dx^ -t- (1 + b)dy^ 4- (1 - b)dz^ , (19) 

where 

h = h{x,t) = A sin [2tt{x - t)] . (20) 

This metric satisfies Einstein's equations only to linear 
order in the wave amplitude A, so if the fully nonlinear 
numerical solution is compared to this approximate solu- 
tion, there will be deviations of order that arise from 
our choice of "analytic solution" rather than from numer- 
ical errors. The amplitude A for the Mexico City tests is 
chosen to be 10~^ so that such deviations in the metric 
components gij are below machine precision. However, 
we still observe 0{A^) deviations in the variables Kij and 
Dkij (which have values of order A), even with A = 10^^, 
because the relative error is well above machine precision. 



A. One-Dimensional Sinusoid 

The sinusoidal waveform chosen in Eq. ([^0]) is only a 
weak test for pseudospectral methods because the Fourier 
basis functions defined in Eqs. (fTT)) and exactly re- 
solve Eq. (PO)) at all times using only three basis func- 
tions; the only truncation errors arc those associated 
with time discretization. Therefore, as a more challeng- 
ing test, in Section llV Bl we repeat the plane wave evo- 
lution using a Gaussian-shaped wave. It is nevertheless 
instructive to evolve the sinusoid and study the resulting 
time discretization errors. Since the dynamics involve no 
change in amplitude, but a change in phase, we expect 
the errors to be primarily phase errors, for reasonably 
small time steps. 
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FIG. 3: Phase Error for 1-D Sinusoidal Linear Wave. 

Phase error at t = 25 crossing times for various time steps, 
and several time-stepping algorithms. These tests were all 
run with three points in the x-direction. The dashed line 
indicates the expected accuracy limit due to roundoff error, 
of. Eq. ((26}. 



This loss of temporal accuracy is particularly relevant 
in efforts to simulate sources for gravitational wave ob- 
servations, as the search for signals involves matching 
expected waveforms against observations. If there is sig- 
nificant error in the phase of the expected waveform, the 
overlap will be poor and detection will be more difficult. 
Although a constant overall scaling error in frequency — 
like the one found in this linear problem — could still re- 
sult in detection, more complex situations would likely 
give rise to more complicated errors. The straightfor- 
ward way to handle this problem is to minimize all time- 
stepping error. 

In Fig. [3] we show the convergence of the phase error in 
the evolution of the sinusoidal linear wave. The solution 
is fully resolved on a 3 x 1 x 1 grid. We keep this grid 
fixed, and decrease the size of the time step. Assuming 
that the only error is some phase error 5<j), the evolved 
gzz will be given by 

= 1 - 10-8sin[27r(a; -t) + . (21) 

At integer multiples of the light crossing time for our 
computational domain, this can be written as 

(722 = 1 — j4 [cos sin (27ra:;) 4- sin cos (27rx)] . (22) 

That is, wc can find the phase error easily from the /c = 1 
sine and cosine components of 1722 (which happen to be 
easily accessible in our code). 

For intermediate time-step sizes, we can see conver- 
gence toward zero phase error with decreasing time step. 
As expected, we observe second-order convergence for It- 
erated Crank-Nicholson stepping, and fourth- and sixth- 
order for the appropriate higher-order Runge-Kutta al- 
gorithms. At very small time-step sizes, a new effect is 
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seen, causing the phase error to increase with decreas- 
ing time step. This effect can be understood as machine 
roundoff error accumulating via a random walk process. 

Suppose we have a variable V{t) that is evolved by 
adding the small changes needed to update its value at 
each time step. Each such operation will introduce a 
fractional error xit) caused by the finite machine preci- 
sion. We assume that the standard tinic-stcp size is At, 
and that there are n intermediate operations in each time 
step. After an evolution through time T, the total error 
added in this way will be 



SV 



iT/At 



(23) 



To avoid tracking each individual error contribution, we 
treat x as a random variable taking values in some range, 
with some probability distribution. 

We estimate the accumulated error due to finite ma- 
chine precision by taking suitable averages over an en- 
semble of random xit) and over a time interval T. If 
there were no asymmetry between positive and negative 
values of x{t), this accumulated error would be zero. Of 
course, we expect almost never to see this case: the most 
likely outcome is an accumulated error comparable to the 
dispersion: 



|<5V| ~ V(5V)2 ^ V(ife)x(t,)x(ife) , (24) 

where the overbar indicates the average over the ensemble 
of random errors xi^)- We can simplify this expression 
by assuming that x{t) has no correlations between time 
steps, and further assuming that the probability distri- 
bution is constant in time and uniform, taking values in 
the range [— e, e], where e is the machine precision. This 
means that x(<j)x(ifc) = Sjk.e^/S. Finally, we approxi- 
mate the discrete time sum as an integral, and obtain 
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^1 



(25) 



We can test this formula by observing its effects in 
the case of phase error for the linear wave. Here, the 
only nontrivial evolved variable is V = gzz, which is 
very nearly 1; so the integral in Eq. (|25p becomes sim- 
ply the evolution time T, which has the value 25 for 
the results plotted in Fig. [H If phase errors dominate, 
Sq^^ AsinSd), so we have 



25 n 
3At 



10- 



(lO-ie) 



10" 



A 



(26) 



where n is assumed to be of order 10. This expression 
is plotted as the dashed hue in Fig. [3l demonstrating 
that the addition of machine-precision errors causes the 
departure from the standard second-, fourth-, and sixth- 
order convergence we observed. From Eq. (j26p we see 
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FIG. 4: Constraints for 1-D Gaussian Linear Wave. 

Here we see the exponential convergence of the constraints 
with higher spatial resolution. At late times, the higher res- 
olutions grow sublinearly in time, probably because of accu- 
mulated machine roundoff error. 



that \6(j)\ is proportional to the ratio e/A; thus \6(j)\ is so 
large in this case because the wave amplitude is so small, 
A = 10-8. 

The phase error is only so clearly visible in these evo- 
lutions because the full solution is described precisely at 
each moment by the first three basis functions. This 
means that discretization error due to spatial differen- 
tiation is essentially at the level of machine precision. 
Indeed, using more than three points actually degrades 
the quality of these one-dimcnsional sinusoid evolutions. 
Power in higher-order basis functions can only be er- 
ror, and hence will necessarily do worse than the low- 
resolution case. We omit plots of the error energy and 
constraints in the higher-resolution cases, as they are very 
nearly the same as those of the more complicated two- 
dimensional evolutions discussed in Section IIVCI 



B. One-Dimensional Gaussian 

As a more challenging test of pseudospectral methods, 
we repeat the one-dimensional linear wave test using a 
periodic Gaussian-shaped wave: 



b ^ A E 6xp 



2w^ 



(27) 



with A = IQ-'*. The summation over j ensures that the 
function is truly periodic at all times. In practice, j need 
only range over a few, depending on the width of the 
Gaussian. The width chosen here is w = 0.05 to ensure 
that features are not too sharp, while still presenting a 
nontrivial challenge to spectral differentiation. 
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FIG. 5: Error Energy for 1-D Gaussian Linear Wave. 

The solid lines show the error energy at 1/2 crossing times, 
with clearly visible exponential convergence at early times. 
The dashed lines show the error energy at integer crossing 
times for the same resolutions. The smallness of the error 
energy at early times demonstrates the low dispersion of the 
numerical method, as explained in the main text. At later 
times, the error is dominated by the quadratic growth ex- 
plained in the text. 



We find behavior comparable to the sinusoidal case, 
although as expected, more collocation points in the x- 
direction are needed to resolve the solution spatially (but 
we still use only a single point in each of the y- and z- 
dimensions). Note the exponential convergence of the 
constraints with increasing resolution in Fig.U) The con- 
straint growth in the highest resolution runs is slower 
than linear in time, and is probably caused by the ac- 
cumulation of errors due to finite machine precision as 
discussed in Sec. IIV Al 

Figure [5] presents the error energy for this run as the 
solid lines. At early times \\5U\\ decreases with resolu- 
tion exponentially to zero, as one would expect. At late 
times, however, \\5U\\ converges toward a parabola. The 
amplitude of this parabola scales in proportion to . In 
the rest of this section, we will first explain a subtlety 
arising when computing followed by a detailed ex- 

planation of why the terms 0{A^) manifest themselves 
in parabolic behavior of 

The comparison of the computed solution with the an- 
alytic solution is performed at the collocation points. By 
virtue of the transformation Eqs. and the errors 
are initially exactly zero at the collocation points. The 
spatial truncation error is nonzero of course, even at the 
initial time; it manifests itself by a deviation of the trun- 
cated series expansion from the analytic solution between 
collocation points. During the evolution, a linear wave 
will simply travel through the computational domain, re- 
turning to the original position after each light-crossing 
time. Since the spectral method has very small disper- 



sion, the evolved shape remains the same. After each 
light-crossing time, therefore, the evolved solution again 
agrees to very high accuracy with the initial analytic so- 
lution at the collocation points. So, comparing the evolu- 
tion with the analytic solution at integer multiples of the 
light-crossing time and at the collocation points will yield 
differences much smaller than spatial truncation error. ^ 
Therefore, a fair comparison that includes the effects of 
spatial truncation error must not be performed at inte- 
ger light-crossing times. These considerations are evident 
from Fig. O where the soHd lines show the "true" \\5U\\ 
observed with 1 /2 light-crossing interval offset, which suf- 
fices because the number of collocation points is always 
odd. The artificially small error energy observed at every 
complete light-crossing interval is shown as dashed lines, 
confirming the excellent low-dispersion property of our 
method. 

At late times the differences between observation at 
full and 1/2 crossing times are swamped by the parabolic 
growth in Similar parabolic deviations of the evo- 

lution from the solution of the linearized equations are 
observed for the other two linear wave evolutions, the 
1-D and 2-D sinusoids (cf. Fig.©. The growth in \\5U\\ 
appears almost entirely due to growth in the k = Q mode 
of diagonal terms in 5gij . Using evolutions of waves with 
different amplitudes and wavelengths, we have verified 
that this growth is proportional to A^t'^/X'^, where A is 
the amplitude and A the wavelength of the disturbance. 
The constant of proportionality depends directly on the 
KST parameter 71 appearing in Eq. ([?]). This parame- 
ter controls the addition of a term ^iNgijC to the ADM 
evolution equation for Kij. The Hamiltonian constraint 
C is roughly constant in time, and varies as A^ / . Since 
the fc = mode of jiNgijC is roughly jiSijC, the fc = 
mode of Ka will grow linearly with time in proportion 
to 7iA^/A^ for each i. That, in turn, will cause small 
quadratic growth in the k = mode of gu . For the more 
well-behaved cases (highest resolutions for the Gaussians; 
all cases for the sinusoids) this model is an excellent fit 
for the observed error energy. 



C. Two-Dimensional Linear Waves 

The linear wave tests above may be modified by ro- 
tating the coordinates by 7r/4 about the z-axis, which 
gives a plane wave propagating along the x-y diagonal. 
By increasing the size of the domain by a factor of \/2 
in each direction, the rotated solution can be made peri- 
odic while maintaining the same wavelength. This con- 
verts the spacetimc from essentially one-dimensional to 
essentially two-dimensional. The purpose of this test 



^ This is also true when comparing at intervals of 1/Af of a crossing 
time if the number of collocation points in the direction of the 
wave's travel is divisible by N. 
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FIG. 6: Constraints for 2-D Linear Waves. The solid 
lines represent the Gaussian wave, while the dashed lines rep- 
resent the sinusoidal wave. The sinusoid is fully resolved spa- 
tially with 3 points. Going to higher resolutions merely intro- 
duces spatial errors in the unnecessary basis functions, which 
leads to an increase in the constraints with resolution. 
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FIG. 7: Error Energy for 2-D Linear Waves. The solid 
lines represent the Gaussian wave, while dashed lines repre- 
sent the sinusoidal wave, both observed at 1/2 crossing times. 
As in the 1-D case (Fig. [S}, both sets of evolutions converge 
to quadratic growth of the error caused by the Hamiltonian 
constraint, explained in the text. 



is to ensure that the symmetry of the one-dimensional 
version does not hide sources of error (although propa- 
gation along a diagonal retains some symmetries). For 
these tests we use a single collocation point in the z- 
dimension, and we vary the (equal) number of collocation 
points in the x- and y-dimensions. We run these tests 
to t = 1000 — ten times longer than is recommended by 
the Apples with Apples collaboration — to better observe 
the stability properties. As shown in Fig. [6l the con- 
straints for the sinusoidal wave increase with increasing 
X- and y-resolution (still using only a single point in the 
^-direction). The constraints for the Gaussian are very 
nearly the same as in the one-dimensional case. Again, 
the A^t^ jy? growth of is visible, as shown in Fig.H 



V. GAUGE WAVE 

The next series of tests involves a simple but time- 
dependent gauge transformation of Minkowski space, in 
the form of a plane wave. The metric used for this Mexico 
City test has the form 



resolve the sinusoidal waveform using only three collo- 
cation points. This is not true for the gauge-wave test, 
because in this case the extrinsic curvature (one of our 
evolved variables) is not a simple sinusoid. Instead, its 
only nonzero component is 



A cos [2tt{x - t)] 
yrTTsin|27r(x^^~i)J 



(30) 



A. One-Dimensional Gauge Wave 

We ran the one-dimensional test described above us- 
ing a single collocation point in each of the y- and 
z-dimensions, and varying the resolution in the x- 
dimension. We find that for both A = 0.01 and A = OA 
our evolution is stable and convergent. Our error energy 
and constraint violations show no signs of instability and 
are strictly better than the filtered two-dimensional evo- 
lutions discussed below. We omit plots for this test be- 
cause the two-dimensional test is more challenging and 
more discriminating. 



ds^ = -(1 + a)dt^ + (1 + a)dx'^ + dy^ + dz^ 



(28) 



A sin [2tt{x - t)] 



(29) 



Two cases are considered: a low amplitude case A = 0.01, 
and a high amplitude case A = 0.1. This is the first test 
for which the nonlinear terms in the equations play an 
important role. 

For the linear plane wave test in Section [TVl we found 
that because wc use a Fourier basis, we were able to fully 



B. Two-Dimensional Gauge Wave 

A simple rotation of coordinates about the z-axis 
makes the wave described by Eqs. (|28p and ([^^ prop- 
agate along the x-y diagonal, as in the case of the linear 
wave. We use an equal number of collocation points in 
the X- and y-dimensions, and a single collocation point 
in z. 
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FIG. 8: Constraints for High-Amplitude 2-D Gauge 
Wave. Dashed lines indicate the unfiltered behavior; solid 
lines indicate the filtered behavior. Note that, despite an 
effective loss of resolution, filtering greatly improves the sta- 
bility of the evolution. 
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FIG. 9: Error Energy for High- Amplitude 2-D Gauge 
Wave. As in Fig. [S] dashed lines are unfiltered, and solid 
lines are filtered. The growth of the filtered error energy is 
exponential in time. For the highest resolution the time step 
was cut in half {dt = dx/80) to reduce time-discretization 
error to the same level as spatial discretization error. The 
dotted line shows the same evolution with time step dt = 
dx/4:0, which is dominated by time discretization error. 



As for the one-dimensional gauge wave test, we run at 
two different amplitudes: A = 0.01 and A = 0.1. For low 
amplitude, A = 0.01, our evolution of the 2-D gauge wave 
is stable and convergent. Again, we omit plots, as our 
results are strictly better than for the more interesting 
high-amplitude case. 



For high amplitude, A = 0.1, we find an exponen- 
tially growing nonconvergent numerical instability, as 
seen in the curves labeled "unfiltered" in Figs. [5] and [S] 
This instability does not appear for the low-amplitude 
case, nor does it appear for either amplitude in the one- 
dimensional gauge wave test. 

The instability appears to be associated with aliasing 
caused by quadratic nonlinearities in the evolution equa- 
tions; this is a well-known phenomenon that often oc- 
curs when applying spectral methods to nonlinear equa- 
tions [Tsj . The basic mechanism for the instability can 
be understood by considering a truncated spectral expan- 
sion for some variable u{x) in terms of N basis functions 



u{x) 



N-l 

k=0 



(31) 



The correct spectral expansion of the expression 
can be obtained by squaring Eq. (|3ip : for most basis 
functions — including the Fourier series of Eq. psp — this 
yields a sum over a total of 2N, and not just N, ba- 
sis functions. But we keep only N basis functions (not 
2A^) in our expansion, so the k > N coefficients of the 
product must be eliminated. Ideally, these k > N coef- 
ficients should be simply discarded and the k < N coef- 
ficients should remain untouched. But it turns out that 
for the pseudospectral method of evaluating nonlinear 
terms (i.e., Fourier transform to obtain values at spa- 
tial collocation points, square these values, then Fourier 
transform back to spectral space) , the power in the extra 
k > N coefficients of the product does not disappear, but 
instead appears as excess power in some of the k < N 
coefficients ( "aliasing" ) . Repeating this process on each 
time step builds up this excess power and produces the 
instability. 

Fortunately, there is a well-known remedy for instabil- 
ities caused by aliasing in nonlinear terms: suppose that 
the upper half (i.e., those with k > N/2) of the coeffi- 
cients Uk in Eq. (PT|) were all zero. Then the spectral 
expansion of u{x)'^ will have zeroes in all its fc > coef- 
ficients, so there is no aliasing, and hence no instability. 
Therefore, we ensure that all coefficients with k > fccut 
are zero by removing those coefficients from the initial 
data and from the right-hand side of the evolution equa- 
tions. It turns out (see, for example, Chapter 11.5 of 
Ref. p^ ) that for a quadratic nonlinearity, it is sufficient 
to filter with fccut = 2A^/3 (and not fccut = N/2) to elimi- 
nate aliasing. As mentioned in Section [TTl the remaining 
number of nonzero coefficients must be odd, which is en- 
sured by reducing fccut by one if necessary. 

The price we pay for stability via this filtering is that 
we must use a factor of 1.5 more spectral coefficients 
(and collocation points) than without filtering in order 
to achieve the same level of spatial discretization error. 
Hence, we use more points for this test than for the pre- 
vious ones: Ni = 15, 21, 27, and 33 points. This leaves 
the effective resolutions at Ni = 9, 13, 17, and 21 points. 
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FIG. 10: Constraints for Shifted Gauge Wave. Solid 
lines indicate A = 0.5, and dashed lines indicate A = 0.1. For 
both amplitudes we filter out the top 1/3 spectral coefficients. 



which are comparable to the resolutions we use in the un- 
filtered case. We see from Figs.|S]and|n]that filtering dra- 
matically reduces the instability. The initial constraint 
violations in these runs, ||C|| « 10~^^, are at the level of 
the finite machine precision, so increasing the resolution 
causes increased — not decreased — constraint violations. 



C. Shifted Gauge Wave 

We also show the results of a new "shifted gauge wave" 
test suggested for addition to the "Apples with Apples" 
suite [ij]. For this test we evolve flat space with the 
usual Minkowski coordinates (t, x, y, z) transformed to 
coordinates [t, x, y, z) via 



y 

z 



X — 

y , 

z . 



47r 
A 

An 



cos [2n{x — t)] , 
■.[2n{x~t)] , 



(32) 

(33) 

(34) 
(35) 



This test includes the effects of a nonvanishing shift vec- 
tor. Wc use the same computational domain and KST 
parameters as in the standard gauge wave tests above. 
The amplitude suggested in Ref. 14| is A = 0.5. We also 
run simulations with A = 0.1. 

At high amplitude, yl = 0.5, we see exponentially grow- 
ing nonconvergent instabilities. Without filtering, the 
code crashes after just a few crossing times. By filtering 
out the top 1/3 spectral coefficients as described above, 
the evolution can be extended as far as t = 60. No other 
choice of filtering seems to improve this further. We also 



FIG. 11: Error Energy for Shifted Gauge Wave. Solid 
lines indicate A — 0.5, while dashed lines indicate A = 0.1. 
The growth in the A — 0.1 runs is roughly linear in time, 
accelerating to quadratic at later times. The dotted line in- 
dicates the standard time step {dt = dx/AO) with 33 points, 
which is dominated by temporal discretization error, while 
the blue dashed curve uses dt = dx/80. 



run the test with an amplitude of A = 0.1. For this am- 
plitude, the evolutions are stable with filtering but un- 
stable without. Figs. [TUl and [TT] show the constraints and 
error energy for these evolutions. The initial constraint 
violations in these runs, ||C|| ~ 10"^'^, are at the level of 
the finite machine precision, so increasing the resolution 
causes increased, rather than decreased, constraint vio- 
lations. The growth in \\SU\\ seen in Fig. [Tljis linear in 
time for t < 100, becoming quadratic at late times. The 
quadratic in time growth is dominated by time-stepping 
error, which tests show is convergent. (Reducing this er- 
ror to the level of spatial truncation error would require 
a prohibitive amount of computing time at the higher 
resolutions.) 



VI. GOWDY SPACETIME 



The Gowdy spacetimes are dynamic cosmological so- 
lutions that present a serious challenge to any numerical 
relativity code. The Gowdy spacetimes are vacuum cos- 
mological models having two spatial Killing fields (planar 
symmetry) that expand from (or, when time-reversed, 
contract toward) a curvature singularity. Two particular 
examples of these spacetimes with relatively simple ana- 
lytical forms were chosen for the Mexico City tests: one 
in which the spacetime expands away from the singular- 
ity; another in which it collapses toward the singularity. 
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FIG. 12: Constraints for Expanding Gowdy Spacetime. 

At early times, the exponential convergence of spectral meth- 
ods is clearly visible. Soon, however, the evolutions are dom- 
inated by constraints growing roughly as e*^^. 



A. Expanding Form 

The metric chosen for the expanding case is 

ds^ = t-'/^e^i^dt^+dz^)+tiePdx^+er''dy^) , (36) 
where 

P{t, z) = Jo(27rt) cos(27rz) , (37) 

\{t,z) = -27rtJo(27rt)Ji(27rt)cos2(27rz) 

^I-kH^ [Jo (27rt) + Jl(2-Kt)\ , (38) 

Ao = A(l,l/8), and J,i is a Bessel function. Asymptot- 
ically, P approaches zero as time increases, and A in- 
creases linearly with time. Because the metric compo- 
nents are singular at i = 0, the Mexico City test begins 
the evolution at i = 1 and proceeds forward in time. 

The time step At required for numerical stability is 
roughly given by the Courant condition At < Ax/w, 
where Ax is the spacing between collocation points and 
V is the coordinate speed of wave propagation, which in 
this case is the coordinate speed of light. For the Gowdy 
metric the coordinate speed of light in the z-direction is 
constant in time, but in the x- and y-directions it varies 
roug hly like i3/4et/2. Therefore, the maximum allowed 
time step decreases in time like t~^/^e~*/^, so for any 
fixed time step, the evolution will soon become numer- 
ically unstable if there is any perturbation in the x- or 
y-directions. This problem can be circumvented by run- 
ning the simulation with just one point in the transverse 
directions, effectively eliminating any perturbation that 
could seed the instability. 




Coordinate Time t 

FIG. 13: Error Energy for Expanding Gowdy Space- 
time. The error energy converges with increased spatial res- 
olution, but ||(5W||2/|1W||2 grows like 6*/^ 



Another difficulty with evolving the expanding Gowdy 
metric is that the metric components and derivatives be- 
come enormous very quickly. By t ~ 725 the numbers 
become larger than 10^^", so the evolution cannot be 
easily handled using standard 64-bit floating-point arith- 
metic. Our evolutions do not actually crash until t = 700; 
unfortunately errors dominate our evolutions long be- 
fore this time, as seen in Fig. 1131 The normalized error 
energy — along with the constraints shown in Fig. 12 — 
grows roughly as e*/^, and accuracy is completely lost in 
these evolutions by t ~ 150. 



B. Collapsing Form 

The time coordinate in the Gowdy metric given 
above can be transformed so that the initial sin- 
gularity is approached only asymptotically in the 
past. The new time coordinate, r, is defined by 
T ee; i In (t/fc), where c = 0.0021195119214607454, and 
k = 9.6707698127640558. The spacetime can be evolved 
backwards indefinitely without reaching the singularity; 
that is, the time step is chosen to be negative. For pur- 
poses of convenience, the evolution is begun at an initial 
time of T = T(io), whereto = 9.8753205829098263, which 
is a zero of Jo(27rt). 

This evolution is far less challenging than the expand- 
ing case. This is because the lapse function is essentially 
an exponential in r, so that the spacetime is becoming 
less dynamical as the simulation progresses and t be- 
comes more negative. The main challenge in this test is 
resolving the spatial features of the solution. For spectral 
methods, the convergence should be exponential with in- 
creasing resolution, which is indeed the behavior shown 
in Figs. [Hand [13 
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FIG. 14: Constraints for Collapsing Gowdy Spacetime. 

Note that the simulation starts at r ~ 9.875, and proceeds 
backwards. 
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FIG. 16: Comparing Constraints for 1-D Gaussian Lin- 
ear Waves. Hamiltonian constraint norms (ragged curves) 
are much smaller than ||C|| for this test, so by themselves are 
not good diagnostics of constraint violations. 
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FIG. 15: Error Energy for Collapsing Gowdy Space- 
time. The simulation starts at r ~ 9.875, and proceeds 
backwards. 



VII. DISCUSSION 

We have applied the full suite of Mexico City tests 0] — 
modified suitably — to a pseudospectral implementation 
of the KST formulation of Einstein's equations. We have 
also implemented the shifted gauge-wave test suggested 
by Babiuc, et al. [l^ . and suggested a number of mi- 
nor changes to the tests that make them better chal- 
lenges for pseudospectral methods. These tests reveal 
that the KST equations with pseudospectral methods 
demonstrate excellent convergence and accuracy, along 
with very good stability in all but a few cases. We have 
derived a fundamental limit Eq. (j25p for the time step 



accuracy possible in a method-of-lines numerical simula- 
tion, and have shown that our implementation is capable 
of quickly achieving that limit in the simple case of a si- 
nusoidal linear wave. We have also shown that the use 
of filtering is very effective in reducing nonlinear aliasing 
instabilities. 

The Mexico City tests provide a basic set of bench- 
marks for evaluating any numerical relativity code: al- 
lowing direct comparisons between different codes that 
use different numerical techniques and different formula- 
tions of the Einstein Equations. However, the tests in 
their present form make too many implicit assumptions 
about the evolution system and the numerical methods. 
Since the creation of the tests, numerical relativity codes 
have become more diverse: using a variety of improved 
numerical techniques (fixed and adaptive mesh refine- 
ment, higher order finite-differencing, multi-block meth- 
ods, spectral methods) and at least two evolution systems 
(generalized harmonic and BSSN) capable of successfully 
evolving binary black hole spacetimes. 

To accommodate the wide range of numerical meth- 
ods and evolution systems now being used, future tests 
need to be formulated in more abstract terms. We rec- 
ommend the following specific changes to the statement 
of the tests: 

1. A code should demonstrate convergence, both spa- 
tial and temporal, appropriate for the numerical 
method used, for each of the tests (gauge wave, lin- 
ear wave, Gowdy spacetime, etc.). 

The number of grid points or the time step needed to 
achieve a given accuracy is highly dependent on the nu- 
merical implementation. Therefore, the test specifica- 
tions should not dictate a certain number of grid points 
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or a certain time-step size as the original formulation of 
the tests did. 

2. The combined error of all evolution variables, and 
the combined constraint violation (including all 
constraints of an evolution system), cf. Eqs. ([7| 
and ([8]), should be reported for each of the tests. 

Prescriptions for examining errors of particular variables 
or constraints, such as those given in the original Mexico 
City tests, are not applicable to evolution systems that do 
not evolve those particular variables or constraints (e.g., 
tetrad or generalized harmonic evolution systems). In 
addition such prescriptions may not encompass all vari- 
ables or constraints (as in the KST system), and may 
therefore fail to detect errors that accumulate only in a 
subset of the evolved variables. To illustrate this point. 
Fig. [TBI shows both the total constraint energy ||C||, and 
the Hamiltonian constraint for the Gaussian linear wave 
(cf. Fig. [4]). The Hamiltonian constraint turns out to be 
anomalously small for the KST system in this case, and 
so is not a good overall error indicator. 

3. Use periodic Gaussian wave spatial profiles in the 
linear and gauge wave tests. 

The sinusoidal spatial profiles specified in the original 
Mexico City tests with periodic boundary conditions pro- 
vide an artificial advantage for spectral techniques. Pe- 
riodic Gaussian profiles arc no more difficult for finite 
difference codes, and provide a significantly greater chal- 
lenge for spectral methods. Finally, 



4. Output data at generic times, not at integer multi- 
ples of the light-crossing time. 

Outputting data at exact integer multiples of the light- 
crossing time significantly underestimates the errors in 
codes with very small dissipation (such as spectral codes) . 

We believe these recommendations will make it easier 
to apply the Mexico City tests fairly to a far wider class 
of numerical relativity codes, and so facilitate apples- 
with-applcs comparisons between these codes. We have 
learned a great deal about the subtle properties of our 
code by carefully running and analyzing these simple 
tests. We encourage other groups to make their results 
from these tests public so that meaningful and objective 
comparisons between codes can be made. 
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